Y=zeros(1,30000);
U=zeros(1,30000);%定义零数组，大小为30000代表30s
wn = 2 ;

ElbsH = [0.2,0.4,0.6,0.8,1.0,2.0];% ξ数组
for L= 1:1:length(ElbsH)%进行循环，对每个ξ值作图
    Elbs = ElbsH(L);
    yk_2 = 0;
    yk_1 = 0;%零初始条件
    uk=0;
    DeltT=0.001;%时间间隔为1ms

    for i=1:1:30000
        if (i>=1000)
            uk=1;%1s时阶跃
        end
        yk=((2+2*Elbs*wn*DeltT)*yk_1 - yk_2 + wn^2 * (DeltT^2) *uk )/(1+ 2*Elbs*wn*DeltT + wn^2 * (DeltT^2));

        Y(i)= yk;%保存yk
        U(i)= uk;%保存uk
        yk_2 = yk_1;
        yk_1 = yk;
    end

    t =  (0:0.001:29.9999);%横坐标转换为时间
    hold on
    plot(t,Y)
end

